body .main-container {
max-width: 100% !important;
width: 100% !important;
}
body {
max-width: 100% !important;
}
body, td {
font-size: 16px;
}
code.r{
font-size: 14px;
}
pre {
font-size: 14px
}knitr::opts_chunk$set(message = FALSE,
warning = FALSE,
fig.align = "center",
out.width = "100%"
)
set.seed(1982)library(INLA)
# Load the data
data(Seeds)
# Create a data frame
df <- data.frame(y = Seeds$r, Ntrials = Seeds$n, Seeds[, 3:5])
# explore the data
df |> head() |> knitr::kable(caption = "Data")| y | Ntrials | x1 | x2 | plate |
|---|---|---|---|---|
| 10 | 39 | 0 | 0 | 1 |
| 23 | 62 | 0 | 0 | 2 |
| 23 | 81 | 0 | 0 | 3 |
| 26 | 51 | 0 | 0 | 4 |
| 17 | 39 | 0 | 0 | 5 |
| 5 | 6 | 0 | 1 | 6 |
The model specification is
\[ y_i\sim\text{Binomial}\left(N_i, p_i\right)\label{eq1}\tag{A} \]
\[ p_i = \text{invlogit}(\eta_i) = \frac{\exp(\eta_i)}{1+\exp(\eta_i)} = \frac{1}{1+\exp(-\eta_i)}\qquad \left(\text{i.e., }\eta_i = \text{logit}(p_i) =\log\left(\frac{p_i}{1-p_i}\right) \right) \]
\[ \eta_i = \beta_0 + \beta_1x_{1i} + \beta_2x_{2i} +v_i \]
\[ \beta_i \sim \text{N}(0, 1000) \]
\[ v_i\sim\text{N}(0,\sigma^2) \] \[ \sigma\sim\text{Exp}(\lambda) \qquad \left(\text{i.e., }\pi(\sigma) = \lambda e^{-\lambda\sigma}\right) \]
\[ \lambda \text{ is such that } \pi(\sigma>1) = 0.01\qquad\left(\text{i.e., }0.01 = \pi(\sigma>1) = 1- \pi(\sigma\leq1) = e^{-\lambda}\implies \lambda = -\dfrac{\log(0.01)}{1}\approx4.6\right)\label{eq2}\tag{B} \]
Expressions from \(\eqref{eq1}\) to \(\eqref{eq2}\) define the model.
list(theta = list(prior = "pc.prec", param = c(1,0.01)))x1x2f()model = "iid" tells inla() we are using
random effectsplate contains the identity of each
observations. That is, plate = 1:nrow(df).The logit() function is not referenced as
"logit" in
control.family = list(control.link = list(model = "logit")).
This "logit" is an option from INLA.
logit <- function(p){log(p/(1-p))}
invlogit <- function(x){exp(x)/(1 + exp(x))}# call to inla
formula1 = y ~ x1 + x2 + f(plate,
model = "iid",
hyper = list(theta = list(prior = "pc.prec",
param = c(1,0.01))))
res1 = inla(formula = formula1,
data = df,
family = "binomial",
Ntrials = Ntrials,
control.family = list(control.link = list(model = "logit")),
control.predictor = list(compute = TRUE, link = 1)) Here we show how to read the summary
#str(res1, 1)
summary(res1)## Time used:
## Pre = 0.375, Running = 0.155, Post = 0.0129, Total = 0.543
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -0.387 0.181 -0.739 -0.390 -0.017 -0.390 0
## x1 -0.360 0.231 -0.838 -0.353 0.076 -0.353 0
## x2 1.033 0.221 0.589 1.034 1.469 1.034 0
##
## Random effects:
## Name Model
## plate IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for plate 19.78 32.61 2.96 10.55 84.94 6.06
##
## Marginal log-Likelihood: -68.44
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
res1$summary.random$plate |> head() |> knitr::kable(format = "markdown", caption = "Summary for the random effects $v_i$. Keep in mind that each of the $v_i$ has a marginal posterior distribution")| ID | mean | sd | 0.025quant | 0.5quant | 0.975quant | mode | kld |
|---|---|---|---|---|---|---|---|
| 1 | -0.3048725 | 0.2799168 | -0.9193189 | -0.2777256 | 0.1690210 | -0.1865579 | 1.8e-06 |
| 2 | -0.0859021 | 0.2322271 | -0.5734249 | -0.0753286 | 0.3594508 | -0.0374485 | 2.0e-07 |
| 3 | -0.3302246 | 0.2477875 | -0.8626290 | -0.3116615 | 0.0969325 | -0.2731563 | 1.2e-06 |
| 4 | 0.2247962 | 0.2451080 | -0.2194151 | 0.2075764 | 0.7478825 | 0.1644008 | 6.0e-07 |
| 5 | 0.0559880 | 0.2459764 | -0.4318958 | 0.0503295 | 0.5620322 | 0.0481711 | 2.0e-07 |
| 6 | 0.1082030 | 0.3233224 | -0.5018691 | 0.0848421 | 0.8176969 | 0.0357540 | 1.6e-06 |
res1$summary.fitted.values |> head() |> knitr::kable(format = "markdown", caption = "Summary for the fitted values $p_i$. Keep in mind that each of the $p_i$ has a marginal posterior distribution")| mean | sd | 0.025quant | 0.5quant | 0.975quant | mode | |
|---|---|---|---|---|---|---|
| fitted.Predictor.01 | 0.3366888 | 0.0585931 | 0.2213399 | 0.3376109 | 0.4489894 | 0.3417845 |
| fitted.Predictor.02 | 0.3851974 | 0.0491712 | 0.2903125 | 0.3845965 | 0.4840581 | 0.3838758 |
| fitted.Predictor.03 | 0.3300258 | 0.0470660 | 0.2394589 | 0.3297444 | 0.4221690 | 0.3291541 |
| fitted.Predictor.04 | 0.4597625 | 0.0579002 | 0.3544738 | 0.4570430 | 0.5794844 | 0.4501589 |
| fitted.Predictor.05 | 0.4190118 | 0.0580471 | 0.3108036 | 0.4165443 | 0.5398323 | 0.4108377 |
| fitted.Predictor.06 | 0.6748795 | 0.0746330 | 0.5208798 | 0.6752715 | 0.8202451 | 0.6717140 |
res1$summary.linear.predictor |> head() |> knitr::kable(format = "markdown", caption = "Summary for the linear predictors $\\eta_i$. Keep in mind that each of the $\\eta_i$ has a marginal posterior distribution")| mean | sd | 0.025quant | 0.5quant | 0.975quant | mode | kld | |
|---|---|---|---|---|---|---|---|
| Predictor.01 | -0.6900710 | 0.2694953 | -1.2578753 | -0.6739590 | -0.2047550 | -0.6301495 | 7e-07 |
| Predictor.02 | -0.4725910 | 0.2100300 | -0.8938666 | -0.4700832 | -0.0637893 | -0.4697858 | 1e-07 |
| Predictor.03 | -0.7160160 | 0.2160388 | -1.1556482 | -0.7093413 | -0.3138757 | -0.7088185 | 3e-07 |
| Predictor.04 | -0.1633026 | 0.2358492 | -0.5994313 | -0.1722528 | 0.3206572 | -0.2093541 | 3e-07 |
| Predictor.05 | -0.3313657 | 0.2414093 | -0.7963651 | -0.3369757 | 0.1596674 | -0.3380490 | 2e-07 |
| Predictor.06 | 0.7522351 | 0.3550642 | 0.0835679 | 0.7321254 | 1.5180090 | 0.6785054 | 1e-07 |
data.frame(p = res1$summary.fitted.values$mean,
eta = res1$summary.linear.predictor$mean,
p.from.eta = invlogit(res1$summary.linear.predictor$mean),
eta.from.p = logit(res1$summary.fitted.values$mean)) |>
head() |>
knitr::kable(format = "markdown", caption = "This shows that `res1$summary.fitted.values` corresponds to the variable ($p_i$) related to the linear predictor ($\\eta_i$) via the link function (`logit()`), and `res1$summary.linear.predictor` corresponds to the linear predictor ($\\eta_i$)")| p | eta | p.from.eta | eta.from.p |
|---|---|---|---|
| 0.3366888 | -0.6900710 | 0.3340173 | -0.6780850 |
| 0.3851974 | -0.4725910 | 0.3840032 | -0.4675454 |
| 0.3300258 | -0.7160160 | 0.3282709 | -0.7080684 |
| 0.4597625 | -0.1633026 | 0.4592648 | -0.1612990 |
| 0.4190118 | -0.3313657 | 0.4179084 | -0.3268315 |
| 0.6748795 | 0.7522351 | 0.6796655 | 0.7303383 |
# betas
res1$summary.fixed |> head() |> knitr::kable(format = "markdown", caption = "Summary for the fixed effects")| mean | sd | 0.025quant | 0.5quant | 0.975quant | mode | kld | |
|---|---|---|---|---|---|---|---|
| (Intercept) | -0.3869998 | 0.1812751 | -0.7387982 | -0.3901860 | -0.0174492 | -0.3900842 | 0e+00 |
| x1 | -0.3598794 | 0.2306549 | -0.8379333 | -0.3525628 | 0.0760315 | -0.3529318 | 1e-07 |
| x2 | 1.0328605 | 0.2214283 | 0.5891781 | 1.0343050 | 1.4690759 | 1.0344034 | 0e+00 |
# hyperparameters
rbind(res1$summary.hyperpar, res1$internal.summary.hyperpar, "exp(Log precision for plate)" = exp(res1$internal.summary.hyperpar)) |>
knitr::kable(format = "markdown", caption = "Summary for the hyperparameters")| mean | sd | 0.025quant | 0.5quant | 0.975quant | mode | |
|---|---|---|---|---|---|---|
| Precision for plate | 19.776422 | 32.6088461 | 2.962334 | 10.552322 | 84.934835 | 6.061489 |
| Log precision for plate | 2.453651 | 0.8231339 | 1.085977 | 2.356346 | 4.441884 | 2.167514 |
| exp(Log precision for plate) | 11.630734 | 2.2776265 | 2.962334 | 10.552322 | 84.934835 | 8.736543 |
This is how the output of inla() corresponds to our
model specification.
\(p_i\) =
res1$summary.fitted.values (just the summary, not the whole
the marginal posterior)
\(\eta_i\) =
res1$summary.linear.predictor (just the summary, not the
whole the marginal posterior)
\(\beta_i\) =
res1$summary.fixed (just the summary, not the whole the
marginal posterior)
\(\tau\) =
res1$summary.hyperpar (just the summary, not the whole the
marginal posterior)
Let us get the marginal posterior for \(\beta_1\). No transformation is needed in this case.
marginal.posterior.beta1.few = res1$marginals.fixed$x1 # we just get a few
marginal.posterior.beta1.all = inla.tmarginal(fun = function(x) x, marginal = res1$marginals.fixed$x1) # we get more
# let us check
data.frame(mean.from.summary = res1$summary.fixed$mean[2],
mean.from.marginal.posterior.few = mean(marginal.posterior.beta1.few[,1]),
mean.from.marginal.posterior.all = mean(marginal.posterior.beta1.all[,1])) |> t() |>
knitr::kable(format = "markdown", caption = "Checking that we are getting the marginal posterior correctly")| mean.from.summary | -0.3598794 |
| mean.from.marginal.posterior.few | -0.3670172 |
| mean.from.marginal.posterior.all | -0.3602327 |
plot(marginal.posterior.beta1.all, type = "p", xlab = expression(beta[1]), ylab = "Probability density", col = "red")
points(marginal.posterior.beta1.few, col = "blue")
legend("topleft", pch = 1, col = c("red", "blue"), legend = c("all samples", "few samples"))To get the marginal posterior in a meaningful scale, we should change
fun from inla.tmarginal(fun, marginal) as
appropriate.
Recall that inla() works with precision and that the
internal parameters are in log() scale. So if we want the
marginal posterior for \(\sigma\) we
need to use a transformation as appropriate.
Recall that \(\tau = \sigma^{-2}\).
So if we let \(\tau\) =
res1$marginals.hyperpar$Precision for plate, then to get
\(\sigma = \tau^{-1/2}\), we need to
use the transformation fun(x) = x^(-1/2).
If we let \(\log(\tau)\) =
res1$internal.marginals.hyperpar$Log precision for plate,
then to get \(\sigma = \tau^{-1/2}\),
we need to use the transformation fun(x) = e^(-x/2).
# posterior of sigma
marginal.posterior.sigma.from.tau = inla.tmarginal(fun = function(x) x^(-0.5), marginal = res1$marginals.hyperpar$`Precision for plate`)
marginal.posterior.sigma.from.log.tau = inla.tmarginal(fun = function(x) exp(-x/2), marginal = res1$internal.marginals.hyperpar$`Log precision for plate`)
plot(marginal.posterior.sigma.from.tau, type = "p", xlab = expression(sigma[iid]), ylab = "Probability density", col = "red")
points(marginal.posterior.sigma.from.log.tau, col = "blue")
legend("topright", pch = 1, col = c("red", "blue"), legend = c("sigma from tau", "sigma from log tau"))So if we let \(\tau\) =
res1$marginals.hyperpar$Precision for plate, then to get
\(\tau\), we need to use the
transformation fun(x) = x.
If we let \(\log(\tau)\) =
res1$internal.marginals.hyperpar$Log precision for plate,
then to get \(\tau\), we need to use
the transformation fun(x) = e^x.
# posterior of tau
marginal.posterior.tau.from.prec = inla.tmarginal(fun = function(x) x, marginal = res1$marginals.hyperpar$`Precision for plate`)
marginal.posterior.tau.from.log.prec = inla.tmarginal(fun = function(x) exp(x), marginal = res1$internal.marginals.hyperpar$`Log precision for plate`)
plot(marginal.posterior.tau.from.prec, type = "p", xlab = expression(tau), ylab = "Probability density", col = "red")
points(marginal.posterior.tau.from.log.prec, col = "blue")
legend("topright", pch = 1, col = c("red", "blue"), legend = c("tau from prec", "tau from log prec"))mean(marginal.posterior.tau.from.prec[,1])## [1] 17.96846
mean(marginal.posterior.tau.from.log.prec[,1])## [1] 17.99411
mean(res1$marginals.hyperpar$`Precision for plate`[,1])## [1] 27.84323
mean(exp(res1$internal.marginals.hyperpar$`Log precision for plate`[,1]))## [1] 27.84323
m.plate.7 = inla.tmarginal(fun = function(x) x, marginal = res1$marginals.random$plate$index.7)
plot(m.plate.7, type="l", xlab = "marginal plate nr 7", ylab = "Probability density")Here are different ways to get \(\eta_7 = \beta_0 + \beta_1x_{17} + \beta_2x_{27} + v_7\)
res1$summary.linear.predictor$mean[7]## [1] 0.8104082
res1$summary.fixed$mean[1] + res1$summary.fixed$mean[2]*df[7,3] + res1$summary.fixed$mean[3]*df[7,4] + mean(res1$marginals.random$plate$index.7[,1])## [1] 0.8140505
res1$summary.fixed$mean[1] + res1$summary.fixed$mean[2]*df[7,3] + res1$summary.fixed$mean[3]*df[7,4] + res1$summary.random$plate$mean[7]## [1] 0.8111264